module poly

use, intrinsic :: iso_c_binding
use basic
implicit none

interface polyval
    module procedure polyvals, polyvalv
end interface ! polyval

contains

subroutine polyvals(p,m,x,y) bind(c)

    real*8, intent(in) :: p(:), x
    integer, intent(in) :: m
    real*8, intent(out) :: y

    integer :: i, n

    n = idnz(p)

    y = 0.D0
    do i = n, m, 1
        y = x*y + p(i)
    end do

    return

end subroutine polyvals

subroutine polyvalv(p,m,x,y,n) bind(c)

    real*8, intent(in) :: p(:), x(:)
    integer, intent(in) :: m, n
    real*8, intent(out) :: y(:)
    integer :: i, j, k

    k = idnz(p)

    !$OMP PARALLEL DEFAULT(NONE), PRIVATE(i,j), SHARED(p,m,x,y,n,k)
    !$OMP DO
    do i = 1, n, 1
        y(i) = 0.d0
        do j = k, m, 1
            y(i) = x(i)*y(i) + p(j)
        end do
    end do
    !$OMP END PARALLEL

    return
end subroutine polyvalv

subroutine polyder(p,m,n,r,l) bind(c)

    real*8, intent(in) :: p(:)
    integer, intent(in) :: m,n
    real*8, intent(inout) :: r(:)
    integer, intent(inout) :: l

    integer :: i, j, k

    if ( l<0 ) then
        l = m - n
        if ( l<=0 ) l = 1
        return
    end if

    if ( n>=m ) then
        r(1) = 0.d0
        return
    end if

    do i = 1, l, 1
        r(i) = p(i)
        k = m-i
        do j = 1, n, 1
            r(i) = r(i)*k
            k = k-1
        end do
    end do
    return

end subroutine polyder

subroutine polyint(p,m,n,r,l) bind(c)

    real*8, intent(in) :: p(:)
    integer, intent(in) :: m,n
    real*8, intent(inout) :: r(:)
    integer, intent(inout) :: l

    integer :: i, j, k

    if ( l<0 ) then
        l = m + n
        return
    end if

    r(m+1:l) = 0.d0

    do i = 1, m, 1
        r(i) = p(i)
        k = m-i
        do j = 1, n, 1
            r(i) = r(i)/dble(k+1)
            k = k+1
        end do
    end do
    return
end subroutine polyint

subroutine polymul(p,m,q,n,r,l) bind(c)

    real*8, intent(in) :: p(:),q(:)
    integer, intent(in) :: m,n
    real*8, intent(inout) :: r(:)
    integer, intent(inout) :: l

    integer :: i, j

    if ( l<0 ) then
        l = m + n - 1
        return
    end if

    do i = 1, l, 1
        r(i) = 0.d0
        do j = 1, m, 1
            if (i-j > -1 .and. i-j < n) &
                r(i) = r(i) + p(i-j+1)*q(j)
        end do
    end do

    return

end subroutine polymul

subroutine polydiv(p,m,q,n,r,l) bind(c)

    real*8, intent(in) :: p(:),q(:)
    integer, intent(in) :: m,n
    real*8, intent(inout) :: r(:)
    integer, intent(inout) :: l

    integer :: i, j

    if ( l<0 ) then
        l = m + n - 1
        return
    end if

    do i = 1, l, 1
        r(i) = 0.d0
        do j = 1, m, 1
            if (i-j > -1 .and. i-j < n) &
                r(i) = r(i) + p(i-j+1)*q(j)
        end do
    end do

    return

end subroutine polydiv

subroutine polyadd(p,m,q,n,r,l) bind(c)

    real*8, intent(in) :: p(:),q(:)
    integer, intent(in) :: m,n
    real*8, intent(inout) :: r(:)
    integer, intent(inout) :: l

    integer :: i, j

    if ( l<0 ) then
        l = max(m,n)
        return
    end if

    if ( m>=n ) then
        do i = 1, l, 1
            if ( i<=m-n ) then
                r(i) = p(i)
            else
                r(i) = p(i) + q(i-m+n)
            end if
        end do
    else
        do i = 1, l, 1
            if ( i<=n-m ) then
                r(i) = q(i)
            else
                r(i) = q(i) + p(i-n+m)
            end if
        end do
    end if

    return

end subroutine polyadd

subroutine polyroots(p,m,r,l) bind(c)
    !> Find roots of the (m-1)-order polynomial by solving the eigen
    ! problem of the coresponding companion matrix.
    !> The roots are all in complex form.

    real*8, intent(in) :: p(:)
    integer, intent(in) :: m
    complex*16, intent(inout) :: r(:)
    integer, intent(inout) :: l

    real*8, allocatable :: A(:,:), work(:), wr(:), wi(:)
    real*8 :: temp(1)
    integer :: i, j, info, lwork

    if ( l<0 ) then
        l = m - 1
        if ( l<1 ) stop "Error: no roots; order too low (<1)!"
        return
    end if

    if ( m == 2 ) then
        r(1) = -p(2)/p(1)
        return
    end if

    allocate(A(l,l))
    allocate(wr(l))
    allocate(wi(l))

    A = 0.d0
    A(1,1) = -p(2)/p(1)
    do i = 2, l, 1
        A(i,i-1) = 1.d0
        A(1,i) = -p(i+1)/p(1)
    end do

    lwork = -1
    call dgeev ( "N", "N", l, A, l, wr, wi, temp, l, temp, l, temp, lwork, info )
    lwork = int(temp(1))
    allocate(work(lwork))
    call dgeev ( "N", "N", l, A, l, wr, wi, temp, l, temp, l, work, lwork, info )

    do i = 1, l, 1
        r(i) = cmplx(wr(i),wi(i))
    end do

    deallocate(A)
    deallocate(work)
    deallocate(wr)
    deallocate(wi)

    return

end subroutine polyroots


function idnz(p) bind(c)
    !> Return the index of first non-zero element in the array
    real*8, intent(in) :: p(:)
    integer :: idnz

    idnz = 1

    do while ( p(idnz)==0.d0 )
        idnz = idnz + 1
    end do
    return
end function idnz

subroutine polyfit(a,t,n,c,oh,ol) bind(c)

    integer, intent(in) :: n
    real(8), intent(in) :: a(:), t(:)
    integer, intent(in) :: oh,ol
    real(8), intent(out) :: c(:)

    integer :: i, j, k
    real(8), allocatable :: b(:)
    real(8), allocatable :: M(:,:)

    k = oh-ol+1
    allocate(M(n,k))
    allocate(b(n))

    c = 0.D0
    b = a

    !$OMP PARALLEL DEFAULT(NONE), PRIVATE(i,j), SHARED(M,t,n,k,oh)
    !$OMP DO
    do i = 1, n, 1
        do j = 1, k, 1
            M(i,j) = t(i)**(oh-j+1)
        end do
    end do
    !$OMP END PARALLEL

    call leastsqs(M,b,n,k)

    c(1:k) = b(1:k)

    deallocate(M)
    deallocate(b)

    return

end subroutine polyfit

subroutine polydetrend(a,n,oh,ol) bind(c)

    integer(4), intent(in) :: n
    real(8), intent(inout) :: a(n)
    integer(4), intent(in) :: oh,ol

    integer :: i
    real(8), allocatable :: t(:), c(:), ac(:)

    if ( oh == 0 ) then
        a = a - sum(a)/dble(n)
        return
    end if

    allocate(t(n))
    allocate(ac(n))
    allocate(c(oh+1))

    t = [ (((i-1)*1.D0),i=1,n) ]
    call polyfit(a,t,n,c,oh,ol)
    call polyval(c,oh+1,t,ac,n)

    a = a - ac

    deallocate(t)
    deallocate(c)
    deallocate(ac)

    return

end subroutine polydetrend

end module poly